##### stimulate 4 * 4 hawkes process with varying baseline intensities

library(data.table)
library(xts)
library(dplyr)
library(pracma)
library(KScorrect)
library(parallel)
library(gmodels)
library(ggplot2)
library(pbmcapply)
setwd("~/Desktop/orderbook/Research/")
source("R/gof/hawkes_utilities.R")

#### Initialize stimulating functions for each single stimulator. i.e. how a arrival of an event stimulates all other events. 
#idx is the stimulator, t is the time after the stimulator event arrives 
phi_mat_func_idx <- function(t, idx){
  if((idx==1) || (idx=="event_1") ){
    ret <- rep(t <=1,6) / 20
  }else if((idx==2) || (idx=="event_2") ){
    ret <- rep(t <=1,6) / 20
  }else if((idx==3) || (idx=="event_3") ){
    ret <- rep(t <=1,6) / 20
#    ret <- rep(0,6) / 10
  }else if((idx==4) || (idx=="event_4") ){
    ret <- rep(t <=2,6) / 20
  }else if((idx==5) || (idx=="event_5") ){
    ret <- rep(t <=2,6) / 20
  }else if((idx==6) || (idx=="event_6") ){
    ret <- rep(t <=2,6) / 20
#    ret <- rep(0,6) / 10
  }
  
  names(ret) <- rep(c("event_1", "event_2", "event_3", "event_4","event_5", "event_6"), each=length(t))
  ret
}

#### a function that combines all stimulators. 
#### Output column is stimulator and row is stimulatee
phi_mat_func <- function(t){
  sapply(1:6, function(idx) phi_mat_func_idx(t, idx))
}


##### initialize all begining parameters 
mu = 0.1/2
tnow = 0
rate_mu = rep(mu, 6)
rate_now = rate_mu
event = data.frame(time=0, event_1=0, event_2=0, event_3 = 0, event_4 = 0, event_5 = 0, event_6 = 0)

#### draw the first event 
first_draw = which(rmultinom(1,1, rep(1/6,6))==1)
if (first_draw == 1){
  event$event_1 =1
}else if (first_draw == 2){
  event$event_2 =1
}else if (first_draw == 3){
  event$event_3 = 1
}else if (first_draw == 4){
  event$event_4 = 1
}else if (first_draw == 5){
  event$event_5 = 1
}else if (first_draw == 6){
  event$event_6 = 1
}

#### rate_sup is the hawkes rates for all events after elevation, 0.1 is the maximum of state increase.
rate_sup = rate_now + phi_mat_func(0) %*% t(data.matrix(event[,-1])) +0.1/2
#### event_idx is the list of current registered event
event_idx <- which(event[1,-1] > 0)


while(tnow < 20000){
  print(tnow)
  time_lapse = rexp(1, rate = sum(rate_sup))
  tnow = tnow + time_lapse
  # 2 is the maximum time that an stimulator could impact the intensity of stimulatee, current is 2.
  relavent_rows = which(event$time > tnow - 2)
  rate_hawkes = sapply(relavent_rows, function(i) phi_mat_func_idx(tnow - event$time[i], event_idx[i]))
  rate_mu_now =rate_mu
  if(floor(tnow) %% 2 == 0){
    rate_mu_now[1] = rate_mu[1] + 0.1/2
    rate_mu_now[2] = rate_mu[2] + 0.1/2
    rate_mu_now[3] = rate_mu[3] + 0.1/2
  }
  if(floor(tnow) %% 3 == 0){
    rate_mu_now[4] = rate_mu[4] + 0.1/2
    rate_mu_now[5] = rate_mu[5] + 0.1/2
    rate_mu_now[6] = rate_mu[6] + 0.1/2
  }
  if(length(relavent_rows) == 0){
    rate_now = rate_mu_now
  }else{
    rate_now = rowSums(rate_hawkes) + rate_mu_now 
  }
  rejection_sample_u = runif(1, 0, sum(rate_sup))
  stopifnot(sum(rate_now) <= sum(rate_sup) + 1e-6)
  event_id = which(rejection_sample_u < cumsum(rate_now))
  if(length(event_id) == 0){
    # 0.1 is value added to make sure rate_sup is strictly larger than rate now at 2/3 elevated points.
    rate_sup = rate_now + 0.1/2
    next
  }
  event_id = event_id[1]
  event_now = data.frame(time=tnow, event_1=0, event_2=0, event_3 =0, event_4 = 0,event_5 =0, event_6 = 0)
  event_now[, 1+event_id] = 1
  event_idx = c(event_idx, event_id)
  event = rbind(event, event_now)
  #  0.1 is value added to make sure rate_sup is strictly larger than rate now at 2/3 elevated points.
  rate_sup = rate_now + phi_mat_func(0) %*% t(data.matrix(event_now[,-1])) + 0.1/2
}



event <- data.table(event)
state <- data.table(cbind((floor(event$time) %% 2 == 0) * 2, (floor(event$time) %% 2 == 0) * 2 ,(floor(event$time) %% 2 == 0) * 2,
                          (floor(event$time) %% 3 == 0) * 3, (floor(event$time) %% 3 == 0) * 3,(floor(event$time) %% 3 == 0) * 3))
colnames(state) <- colnames(event)[-1]

state_func <- function(t){
  cbind((floor(t) %% 2 == 0) * 2, (floor(t) %% 2 == 0) * 2, (floor(t) %% 2 == 0) * 2,
        (floor(t) %% 3 == 0)*3, (floor(t) %% 3 == 0)*3, (floor(t) %% 3 == 0)*3)
}

delta = 0.1
support_max = 4
support_seq = seq(0,support_max,delta)

hawkes <- estimate_hawkes(event = data.table(event), delta = delta, support_max = support_max, state_func = state_func)
hawkes$phi_mat_func <- phi_mat_func_idx
jpeg(file = "~/Desktop/orderbook/Results/simulation_6*6.png")
plot_hawkes(hawkes,matrix_layout=TRUE)
dev.off()
head(hawkes$hawkes_kernel_est)


par(mfrow=c(6,6))
par(mar=c(1.1,1,1,1))

time_scale <- (1:floor(hawkes$support_max/hawkes$delta)) * hawkes$delta
for(stimulater in colnames(hawkes$hawkes_kernel_est)){
  for(stimulatee in colnames(hawkes$hawkes_kernel_est)){
    kernel_this <- hawkes$hawkes_kernel_est[rownames(hawkes$hawkes_kernel_est)==paste0("design_matrix", stimulater),]
    if(!is.null(hawkes$phi_mat_func)){
      phi_func <- function(t) {
        ret <- phi_mat_func_idx(t, stimulater)
        ret[names(ret) == stimulatee]
      }
    }else{
      phi_func <- NULL
    }
    plot_hawkes_kernel(rev(time_scale), kernel_this[,stimulatee], phi_func)
    title(paste0("event of ",stimulater," stimulate ",stimulatee))
    
  }
}








##### Generate QQ-plot ##########

event = data.table(event)
event_list = list()
for (i in 1:length(event$event_1)){
  if (event$event_1[i] == 1){
    event_list[[i]] = "event_1"
  }else if(event$event_2[i] == 1){
    event_list[[i]] = "event_2"
  }else if(event$event_3[i] == 1){
    event_list[[i]] = "event_3"
  }else if(event$event_4[i] == 1){
    event_list[[i]] = "event_4"
  }else if(event$event_5[i] == 1){
    event_list[[i]] = "event_5"
  }else if(event$event_6[i] == 1){
    event_list[[i]] = "event_6"
  }
}
event[,event_type:= unlist(event_list)]

event_type = c("event_1", "event_2", "event_3", "event_4", "event_5", "event_6")
est = hawkes$hawkes_kernel_est
est[is.na(est)] = 0
est_events = est[rownames(est) %like% "design",]
est_state = est[rownames(est) %like% "state",]


get_phi <- function(interval, stimulator, stimulatee){
  if(near(interval,support_max)==TRUE){
    phi = 0
  }
  else{
    float_epsilon = near(interval, support_seq)
    if(sum(float_epsilon) > 0){
      interval = support_seq[which(near(interval, support_seq)==1)]
    }
    index = support_max/delta - floor(interval/delta)
    phi = est_events[rownames(est_events)==paste0("design_matrix", stimulator),stimulatee][[index]]
  }
  phi
}


get_lambda<-function(t, idx=selected_type){
  print (t)
  sub = event%>%filter(between(time, t-support_max,t))
  state_value = state_func(t)
  colnames(state_value) = event_type
  state_value = toString(state_value[1,idx])
  
  if(dim(sub)[1] == 0){
    lambda = est_state[rownames(est_state) %like% state_value ,][[idx]]
  }
  else{
    sub = data.table(sub)
    sub[, ("interval") := t - time]
    sub[,excites := mapply(get_phi, interval, event_type,idx)]
    lambda = sum(sub[,excites]) + est_state[rownames(est_state) %like% state_value ,][[idx]] 
  }
  lambda 
}
##### Generate change time of step function lambda(t)-----------------------------------
change_time = c()
change_time = lapply(1:dim(event)[1], function(i){
  c(change_time, seq(event$time[i], event$time[i]+support_max, delta  ))
})
change_time = unlist(change_time)
state_time_change = 0: floor(max(event$time))
change_time = sort(unique(c(change_time, state_time_change)))

#lambda_list = pbmclapply(change_time, get_lambda, mc.cores = 10)
selected_type = event_type[1]
lambda_list = pbmclapply(change_time, get_lambda, mc.cores = 10)
lambda_list = unlist(lambda_list)



par(mfrow=c(3,2))
for (i in c(1,2,3,4,5,6)){
  print(selected_type)
  selected_type = event_type[i]
  lambda_list = pbmclapply(change_time, get_lambda, mc.cores = 10)
  lambda_list = unlist(lambda_list)
  selected_event = event[event[[selected_type]]==1,1]
  integral = pbmclapply(1:length(selected_event$time), function(i){
    curr_time = selected_event$time[i]
    index = which(near(change_time, curr_time) == TRUE)
    sum(lambda_list[1:index-1] * diff(change_time[1:index]))
  },mc.cores = 10)
  integral = unlist(integral)
  empirical = diff(integral)
  
  exp = rexp(length(empirical))
  
  if(selected_type=="event_1"){
    modified_type = "+1(i)"
  }else if(selected_type == "event_2"){
    modified_type = "+1(c)"
  }else if(selected_type == "event_3"){
    modified_type = "+1(t)"
  }else if(selected_type == "event_4"){
    modified_type = "-1(i)"
  }else if(selected_type == "event_5"){
    modified_type = "-1(c)"
  }else if(selected_type == "event_6"){
    modified_type = "-1(t)"
  }
  
  qqplot(x=exp, y = empirical, main = paste0("QQ Plot for ", modified_type , ", T = 20000s"),
         xlab = "Standard Exponential Distribution",
         ylab = "rescaled Time")
  qqline(empirical, distribution=qexp)
  
}






#selected_type = event_type[3]
selected_type = event_type[2]
lambda_list = pbmclapply(change_time, get_lambda, mc.cores = 10)
lambda_list = unlist(lambda_list)

selected_event = event[event[[selected_type]]==1,1]
integral = pbmclapply(1:length(selected_event$time), function(i){
  curr_time = selected_event$time[i]
  index = which(near(change_time, curr_time) == TRUE)
  sum(lambda_list[1:index-1] * diff(change_time[1:index]))
},mc.cores = 10)
integral = unlist(integral)
empirical = diff(integral)

exp = rexp(length(empirical))
qqplot(x=exp, y = empirical, main = paste0("QQ Plot for ", selected_type , ", T = 20000s"),
       xlab = "Standard Exponential Distribution",
       ylab = "rescaled Time")
qqline(empirical, distribution=qexp)
ks.test(empirical, "pexp",1)





######### Write a function that simulate hawkes process

simulate_hawkes_with_mu <- function(){
  mu = 0.1/2
  tnow = 0
  rate_mu = rep(mu, 6)
  rate_now = rate_mu
  event = data.frame(time=0, event_1=0, event_2=0, event_3 = 0, event_4 = 0, event_5 = 0, event_6 = 0)
  
  #### draw the first event 
  first_draw = which(rmultinom(1,1, rep(1/6,6))==1)
  if (first_draw == 1){
    event$event_1 =1
  }else if (first_draw == 2){
    event$event_2 =1
  }else if (first_draw == 3){
    event$event_3 = 1
  }else if (first_draw == 4){
    event$event_4 = 1
  }else if (first_draw == 5){
    event$event_5 = 1
  }else if (first_draw == 6){
    event$event_6 = 1
  }
  
  #### rate_sup is the hawkes rates for all events after elevation, 0.1 is the maximum of state increase.
  rate_sup = rate_now + phi_mat_func(0) %*% t(data.matrix(event[,-1])) +0.1/2
  #### event_idx is the list of current registered event
  event_idx <- which(event[1,-1] > 0)
  
  
  while(tnow < 20000){
#    print(tnow)
    time_lapse = rexp(1, rate = sum(rate_sup))
    tnow = tnow + time_lapse
    # 2 is the maximum time that an stimulator could impact the intensity of stimulatee, current is 2.
    relavent_rows = which(event$time > tnow - 2)
    rate_hawkes = sapply(relavent_rows, function(i) phi_mat_func_idx(tnow - event$time[i], event_idx[i]))
    rate_mu_now =rate_mu
    if(floor(tnow) %% 2 == 0){
      rate_mu_now[1] = rate_mu[1] + 0.1/2
      rate_mu_now[2] = rate_mu[2] + 0.1/2
      rate_mu_now[3] = rate_mu[3] + 0.1/2
    }
    if(floor(tnow) %% 3 == 0){
      rate_mu_now[4] = rate_mu[4] + 0.1/2
      rate_mu_now[5] = rate_mu[5] + 0.1/2
      rate_mu_now[6] = rate_mu[6] + 0.1/2
    }
    if(length(relavent_rows) == 0){
      rate_now = rate_mu_now
    }else{
      rate_now = rowSums(rate_hawkes) + rate_mu_now 
    }
    rejection_sample_u = runif(1, 0, sum(rate_sup))
    stopifnot(sum(rate_now) <= sum(rate_sup) + 1e-6)
    event_id = which(rejection_sample_u < cumsum(rate_now))
    if(length(event_id) == 0){
      # 0.1 is value added to make sure rate_sup is strictly larger than rate now at 2/3 elevated points.
      rate_sup = rate_now + 0.1/2
      next
    }
    event_id = event_id[1]
    event_now = data.frame(time=tnow, event_1=0, event_2=0, event_3 =0, event_4 = 0,event_5 =0, event_6 = 0)
    event_now[, 1+event_id] = 1
    event_idx = c(event_idx, event_id)
    event = rbind(event, event_now)
    #  0.1 is value added to make sure rate_sup is strictly larger than rate now at 2/3 elevated points.
    rate_sup = rate_now + phi_mat_func(0) %*% t(data.matrix(event_now[,-1])) + 0.1/2
  }
  
  
  
  event <- data.table(event)
  state <- data.table(cbind((floor(event$time) %% 2 == 0) * 2, (floor(event$time) %% 2 == 0) * 2 ,(floor(event$time) %% 2 == 0) * 2,
                            (floor(event$time) %% 3 == 0) * 3, (floor(event$time) %% 3 == 0) * 3,(floor(event$time) %% 3 == 0) * 3))
  colnames(state) <- colnames(event)[-1]
  
  state_func <- function(t){
    cbind((floor(t) %% 2 == 0) * 2, (floor(t) %% 2 == 0) * 2, (floor(t) %% 2 == 0) * 2,
          (floor(t) %% 3 == 0)*3, (floor(t) %% 3 == 0)*3, (floor(t) %% 3 == 0)*3)
  }
  
  delta = 0.1
  support_max = 4
  support_seq = seq(0,support_max,delta)
  
  hawkes <- estimate_hawkes(event = data.table(event), delta = delta, support_max = support_max, state_func = state_func)
  hawkes$phi_mat_func <- phi_mat_func_idx
  return(list(hawkes))
}


hawkes_ret = list()
hawkes_kernel_ret = list()
state_df = data.frame()
for (i in 1:100){
  print(i)
  ret = simulate_hawkes_with_mu()
  hawkes = ret[[1]]
  hawkes_ret[[i]] = hawkes
  hawkes_kernel_ret[[i]] = hawkes$hawkes_kernel_est
  state_df[i,1] = hawkes$hawkes_kernel_est[1,1]
  state_df[i,2] = hawkes$hawkes_kernel_est[1,2]
  state_df[i,3] = hawkes$hawkes_kernel_est[1,3]
  state_df[i,4] = hawkes$hawkes_kernel_est[1,4]
  state_df[i,5] = hawkes$hawkes_kernel_est[1,5]
  state_df[i,6] = hawkes$hawkes_kernel_est[1,6]
  state_df[i,7] = hawkes$hawkes_kernel_est[2,1]
  state_df[i,8] = hawkes$hawkes_kernel_est[2,2]
  state_df[i,9] = hawkes$hawkes_kernel_est[2,3]
  state_df[i,10] = hawkes$hawkes_kernel_est[3,4]
  state_df[i,11] = hawkes$hawkes_kernel_est[3,5]
  state_df[i,12] = hawkes$hawkes_kernel_est[3,6]
#  state_df[i,1] = hawkes$hawkes_kernel_est[1,1]
#  state_df[i,2] = hawkes$hawkes_kernel_est[1,2]
#  state_df[i,3] = hawkes$hawkes_kernel_est[1,3]
#  state_df[i,4] = hawkes$hawkes_kernel_est[1,4]
#  state_df[i,5] = hawkes$hawkes_kernel_est[2,1]
#  state_df[i,6] = hawkes$hawkes_kernel_est[2,2]
#  state_df[i,7] = hawkes$hawkes_kernel_est[3,3]
#  state_df[i,8] = hawkes$hawkes_kernel_est[3,4]
  #  state_df[i,1] = hawkes$hawkes_kernel_est[1,1]
  #  state_df[i,2] = hawkes$hawkes_kernel_est[1,2]
  #  state_df[i,3] = hawkes$hawkes_kernel_est[2,1]
  #  state_df[i,4] = hawkes$hawkes_kernel_est[3,2]
}

colnames(state_df) = c("base_1", "base_2", "base_3", "base_4","base_5", "base_6",
                       "state21", "state22","state23", "state34", "state35","state36")

state_df1 = as.data.frame(t(apply(as.matrix(state_df), 2, function(x) ci(x))))
names(state_df1) = c("mean", "ci_lower", "ci_upper", "sd")
#x_names = c("Event 1 baseline", "Event 2 baseline", "Event 3 baseline", "Event 4 baseline","Event 5 baseline","Event 6 baseline",
#            "Event 1 state + baseline", "Event 2 state + baseline","Event 3 state + baseline",
#            "Event 4 state + baseline", "Event 5 state + baseline","Event 6 state + baseline")

x_names = c("+1(i) baseline", "+1(c) baseline", "+1(t) baseline", "-1(i) baseline","-1(c) baseline","-1(t) baseline",
            "+1(i) state + baseline", "+1(c) state + baseline","+1(t) state + baseline",
            "-1(i) state + baseline", "-1(c) state + baseline","-1(t) state + baseline")



ggplot(state_df1) +
  geom_bar( aes(x=x_names , y=mean), stat="identity", fill="skyblue", alpha=0.7) +
  geom_point(aes(x = x_names, y = c(rep(0.05,6), rep(0.10,6)), colour = "True Value"),
             stat="identity",
             alpha=.8,
             size=3) + 
  geom_errorbar( aes(x=x_names, ymin=ci_lower, ymax=ci_upper), width=0.4, colour="orange", alpha=0.9, size=1) + 
  ylab("Mean estimation") + xlab("")  + labs(color='') +  scale_colour_manual(name="",  
                                                                              values = c("True Value"="black")) +
  scale_x_discrete(guide = guide_axis(angle = 90))

phi_mat_func_idx <- function(t, idx){
  if((idx==1) || (idx=="+1(i)") ){
    ret <- rep(t <=1,6) / 20
  }else if((idx==2) || (idx=="+1(c)") ){
    ret <- rep(t <=1,6) / 20
  }else if((idx==3) || (idx=="+1(t)") ){
    ret <- rep(t <=1,6) / 20
    #    ret <- rep(0,6) / 10
  }else if((idx==4) || (idx=="-1(i)") ){
    ret <- rep(t <=2,6) / 20
  }else if((idx==5) || (idx=="-1(c)") ){
    ret <- rep(t <=2,6) / 20
  }else if((idx==6) || (idx=="-1(t)") ){
    ret <- rep(t <=2,6) / 20
    #    ret <- rep(0,6) / 10
  }
  
  names(ret) <- rep(c("+1(i)", "+1(c)", "+1(t)", "-1(i)", "-1(c)", "-1(t)"), each=length(t))
  ret
}




plot_hawkes_new <- function(hawkes){
#  par(mfrow=c(6,6))
#  par(mar=c(1.1,1,1,1))
  layout(matrix(1:36, byrow=TRUE, ncol=6,nrow=6))
  
  time_scale <- (1:floor(hawkes$support_max/hawkes$delta)) * hawkes$delta
  for(stimulater in colnames(hawkes$hawkes_kernel_est)){
    for(stimulatee in colnames(hawkes$hawkes_kernel_est)){
      kernel_this <- hawkes$hawkes_kernel_est[rownames(hawkes$hawkes_kernel_est)==paste0("design_matrix", stimulater),]
      if(!is.null(hawkes$phi_mat_func)){
        phi_func <- function(t) {
          ret <- phi_mat_func_idx(t, stimulater)
          ret[names(ret) == stimulatee]
        }
      }else{
        phi_func <- NULL
      }
      par(mar=c(2,4,2,2))
      plot_hawkes_kernel(rev(time_scale), kernel_this[,stimulatee], phi_func)
      title(paste0("",stimulater," stimulate ",stimulatee))
      
    }
  }
  
}

plot_hawkes_kernel <- function(time_pints, phi_value, phi_true_func=NULL){
  plot(time_pints, phi_value, xlab="", ylab="Hawkes kernel", ylim=range(c(0, phi_value), na.rm = TRUE))
  finite_idx <- which(is.finite(phi_value))
  if(length(finite_idx) >= 4){
    fitted_curve <- smooth.spline(time_pints[finite_idx], phi_value[finite_idx])
    lines(fitted_curve, col=2)
  }
  if(!is.null(phi_true_func)){
    time_true <- seq(0, max(time_pints), length=101)
    lines(time_true, sapply(time_true, phi_true_func), col=4)
  }
}




hawkes_final = hawkes
hawkes_final$hawkes_kernel_est = Reduce("+", hawkes_kernel_ret) / length(hawkes_kernel_ret)
colnames(hawkes_final$hawkes_kernel_est) = c("+1(i)", "+1(c)", "+1(t)", "-1(i)", "-1(c)", "-1(t)")
rownames(hawkes_final$hawkes_kernel_est) = c(rownames(hawkes_final$hawkes_kernel_est)[1:3], 
                                             rep(c("design_matrix+1(i)", "design_matrix+1(c)", "design_matrix+1(t)", 
                                                   "design_matrix-1(i)", "design_matrix-1(c)", "design_matrix-1(t)"), support_max/delta))

plot_hawkes_new(hawkes_final)
head(hawkes_final$hawkes_kernel_est)

saveRDS(hawkes_final, "~/Desktop/orderbook/Research/R/gof/hawkes_final.rds")
saveRDS(hawkes_ret, "~/Desktop/orderbook/Research/R/gof/hawkes_ret.rds")

hawkes_final <- readRDS("~/Desktop/orderbook/Research/R/gof/hawkes_final.rds")
hawkes_ret <- readRDS("~/Desktop/orderbook/Research/R/gof/hawkes_ret.rds")